sRMS <- function(a, ref, samp, time, B, D, lambda, kappa)
{
  w.tmp <- time + t(B %*% a )
  P.tmp <- D %*% a
  P <- lambda * t(P.tmp) %*% P.tmp + kappa * a %*% a
  
  interp <- t(apply(samp, 1, function(x) interpol(w.tmp, x)))
  
  if (nrow(ref) == 1) ref <- c(ref)
  r <- interp - ref

  sqrt(mean(r^2, na.rm = TRUE) + P / length(ref))
}
